
import numpy as np
import xarray as xr


#ds = xr.open_dataset("https://thredds-su.ipsl.fr/thredds/dodsC/ipsl_thredds/brocksce/PALEORCA2.3/14Ma/IPSL.PALEORCA2_gridCF.nc")
ds = xr.open_dataset("https://thredds-su.ipsl.fr/thredds/dodsC/ipsl_thredds/brocksce/PALEORCA2.3/17Ma/IPSL.PALEORCA2_gridCF.nc")
ds

<xarray.Dataset>
Dimensions:       (y: 149, x: 182, four: 4, depth: 31, two: 2)
Coordinates:
    lon           (y, x) float32 ...
    lat           (y, x) float32 ...
  * depth         (depth) float32 5.0 15.0 25.0 ... 4.25e+03 4.75e+03 5.25e+03
Dimensions without coordinates: y, x, four, two
Data variables:
    bounds_lon    (y, x, four) float32 ...
    bounds_lat    (y, x, four) float32 ...
    bounds_depth  (depth, two) float32 ...
    MASK          (depth, y, x) float32 ...
    MASKDRAW      (depth, y, x) float32 ...
    VOL           (depth, y, x) float32 ...
    AREA          (y, x) float32 ...
    BATHY         (y, x) float32 ...
Attributes:
    Conventions:  CF-1.1
    institution:  PALEORCA2
    production:   NEMO

[27118 values with dtype=float32]

[27118 values with dtype=float32]

array([4.999938e+00, 1.500029e+01, 2.500176e+01, 3.500541e+01, 4.501332e+01,
       5.502950e+01, 6.506181e+01, 7.512551e+01, 8.525037e+01, 9.549429e+01,
       1.059699e+02, 1.168962e+02, 1.286979e+02, 1.421953e+02, 1.589606e+02,
       1.819628e+02, 2.166479e+02, 2.724767e+02, 3.643030e+02, 5.115348e+02,
       7.322009e+02, 1.033217e+03, 1.405698e+03, 1.830885e+03, 2.289768e+03,
       2.768242e+03, 3.257479e+03, 3.752442e+03, 4.250401e+03, 4.749913e+03,
       5.250227e+03], dtype=float32)

[108472 values with dtype=float32]

[108472 values with dtype=float32]

[62 values with dtype=float32]

[840658 values with dtype=float32]

[840658 values with dtype=float32]

[840658 values with dtype=float32]


# MASKDRAW or MASK can be used since duplicate polys will be removed later
#mask = xr.where(ds['MASKDRAW'].isel(depth=0) == 1, 1, np.nan)
mask = xr.where(ds['MASK'].isel(depth=0) == 1, 1, np.nan)
mask.plot()

<matplotlib.collections.QuadMesh at 0x7f0bea10a3b0>


# build boundaries and and mask them

mask = xr.where(ds['MASK'].isel(depth=0) == 1, 1, np.nan)

mask = mask.to_numpy().reshape(-1)
print(mask.shape)

mask1 = mask[~np.isnan(mask)]
#mask1 = mask
print(mask1.shape[0])

blon = ds['bounds_lon'].to_numpy()
blat = ds['bounds_lat'].to_numpy()
nvertex = blon.shape[-1]
print(nvertex)

blon = blon.reshape(-1, nvertex)
blat = blat.reshape(-1, nvertex)
print(blon.shape)

mask = np.repeat(mask, nvertex).reshape(mask.shape[0], nvertex)

blon = blon[~np.isnan(mask)].reshape(mask1.shape[0], nvertex)
blat = blat[~np.isnan(mask)].reshape(mask1.shape[0], nvertex)
print(blon.shape)

(27118,)
16385
4
(27118, 4)
(16385, 4)


deg2rad = np.pi/180.
x = np.cos(blat*deg2rad)*np.cos(blon*deg2rad)
y = np.cos(blat*deg2rad)*np.sin(blon*deg2rad)
z = np.sin(blat*deg2rad)


# Structures for polydata
points = np.stack((x,y,z), axis=2).reshape(x.size, 3)
faces = np.arange(x.shape[0] * nvertex).reshape(x.shape[0], nvertex)
faces = np.insert(faces, 0, nvertex, axis=1)


import pyvista as pv


mesh = pv.PolyData(points, faces)
mesh


# Clean polydata with tolerance set
mesh = pv.PolyData(points, faces)
mesh.clean(inplace=True, tolerance=0.001)
mesh


import vtk


# Duplicate cells removed
removePolyData = vtk.vtkRemoveDuplicatePolys()
removePolyData.SetInputData(mesh)
removePolyData.Update()

mesh = removePolyData.GetOutput()
mesh = pv.wrap(mesh)
mesh


edges = mesh.extract_feature_edges(boundary_edges=True)
edges


import pyvista as pv
import vtk

vtk.vtkLogger.SetStderrVerbosity(vtk.vtkLogger.VERBOSITY_OFF)

edges = mesh.extract_feature_edges(boundary_edges=True)

pl = pv.Plotter()
pl.set_background('white')

actor1 = pl.add_mesh(pv.Sphere(radius=0.999, theta_resolution=360, phi_resolution=180))
actor2 = pl.add_mesh(mesh, show_edges=True, edge_color="gray")
actor3 = pl.add_mesh(edges, color="red", line_width=2)
viewer = pl.show(jupyter_backend='panel', return_viewer=True)
display(viewer)


edges.GetNumberOfLines()

1624


edges.connectivity()


regions = edges.connectivity()
regCount = len(set(pv.get_array(regions ,name="RegionId")))
regCount

14


import random


pl = pv.Plotter()
pl.set_background('white')
pl.add_mesh(pv.Sphere(radius=0.999, theta_resolution=360, phi_resolution=180))

CONT_LON = []
CONT_LAT = []

regions = edges.connectivity()
regCount = len(set(pv.get_array(regions, name="RegionId")))

connectivityFilter = vtk.vtkPolyDataConnectivityFilter()
stripper = vtk.vtkStripper()

for r in range(regCount):
    connectivityFilter.SetInputData(edges)
    connectivityFilter.SetExtractionModeToSpecifiedRegions()
    connectivityFilter.InitializeSpecifiedRegionList()
    connectivityFilter.AddSpecifiedRegion(r)
    connectivityFilter.Update()
    
    reg = connectivityFilter.GetOutput()
    
    stripper.SetInputData(reg)
    stripper.SetJoinContiguousSegments(True)
    stripper.Update()
    reg = stripper.GetOutput()
      
    pts = reg.GetCell(0).GetPoints()
    thetas = []
    phis = []
    for i in range(pts.GetNumberOfPoints()):
        x,y,z = pts.GetPoint(i)
        theta   =  90 - np.arccos(z)*180 / np.pi
        phi     =  np.arctan2(y,x)*180 / np.pi
        phis.append(phi)
        thetas.append(theta)
    # as to be reversed https://github.com/SciTools/cartopy/issues/2103
    CONT_LON.extend(reversed(phis))
    CONT_LAT.extend(reversed(thetas))
    CONT_LON.append(-1E34)
    CONT_LAT.append(-1E34)
    
    random_color = "#"+''.join([random.choice('0123456789ABCDEF') for i in range(6)])
    pl.add_mesh(reg, color=random_color, line_width=2)
    
viewer = pl.show(jupyter_backend='panel', return_viewer=True)
display(viewer)


import netCDF4

rootgrp = netCDF4.Dataset("paleo_continents.nc", "w", format="NETCDF4", clobber=True)
npoints = rootgrp.createDimension("npoints", None)
longitudes = rootgrp.createVariable("CONT_LON","f8",("npoints",))
latitudes = rootgrp.createVariable("CONT_LAT","f8",("npoints",))
longitudes[:] = CONT_LON
latitudes[:] = CONT_LAT
longitudes.missing_value = -1.E+34
latitudes.missing_value = -1.E+34
rootgrp.close()


import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import xarray as xr

ds = xr.open_dataset("paleo_continents.nc")

x = ds['CONT_LON'] ; y = ds['CONT_LAT']

fig = plt.figure(figsize=(10,5))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.Robinson(10))

#transform = ccrs.PlateCarree()
transform = ccrs.Geodetic()

#ax.plot(x, y, transform=transform, color='red', lw=0.5)
ax.fill(x, y, transform=transform, color='lightgray', lw=0.5, edgecolor="black")
ax.set_global()
ax.gridlines()
plt.show()

ds.close()


import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import xarray as xr

ds = xr.open_dataset("paleo_continents.nc")

x = ds['CONT_LON'] ; y = ds['CONT_LAT']

fig = plt.figure(figsize=(10,5))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.Robinson(10))

transform = ccrs.Geodetic()

indices = np.where(np.isnan(x.values))[0]
#print(indices)

indice_curr = 0
for indice in indices:
    xselect = x[indice_curr:indice]
    yselect = y[indice_curr:indice]
    random_color = "#"+''.join([random.choice('0123456789ABCDEF') for i in range(6)])
    #ax.plot(xselect, yselect, transform=transform, color=random_color, lw=0.5)
    ax.fill(xselect, yselect, transform=transform, color=random_color, lw=0.5, edgecolor="black")
    indice_curr = indice+1
    

ax.set_global()
ax.gridlines()
plt.show()

ds.close()

PolyData	Information
N Cells	16385
N Points	65540
N Strips	0
X Bounds	-1.000e+00, 1.000e+00
Y Bounds	-1.000e+00, 1.000e+00
Z Bounds	-9.748e-01, 1.000e+00
N Arrays	0

PolyData	Information
N Cells	16385
N Points	17032
N Strips	0
X Bounds	-1.000e+00, 1.000e+00
Y Bounds	-1.000e+00, 1.000e+00
Z Bounds	-9.748e-01, 1.000e+00
N Arrays	0

PolyData	Information
N Cells	16233
N Points	17032
N Strips	0
X Bounds	-1.000e+00, 1.000e+00
Y Bounds	-1.000e+00, 1.000e+00
Z Bounds	-9.748e-01, 1.000e+00
N Arrays	0

PolyData	Information
N Cells	1624
N Points	1623
N Strips	0
X Bounds	-7.855e-01, 9.972e-01
Y Bounds	-9.779e-01, 9.900e-01
Z Bounds	-9.748e-01, 9.941e-01
N Arrays	0

Extracting continent shorelines¶

Example with a paleo orca 2.3 configuration¶